#### Load the data ####

# Impute SE from replication data (temporary)

library(rdrobust)
#library(sandwich)

#### Load data and run analysis for studies that we impute via replication data/code ####

## three articles/estimates require this:
#caughey_etal_2017: compute standard error and p-value
#carson_sievert_2017: compute point estimate, standard error and p-value
#dbk_warshaw_2016: compute standard error

#(1) caughey_etal_2017:

caughey_etal_2017_df <- read.csv("data/cl-replication-data/caughey_etal_2017.csv")

caughey_etal_2017_res <- rdrobust(
  y = caughey_etal_2017_df$GovLibD1, 
  x = caughey_etal_2017_df$DemMarginGov,
  bwselect = "mserd", 
  all = TRUE, 
  kernel = "tri"
)

rep_res <- data.frame(
  "estimate_code" = "caughey_etal_2017",
  "beta_est" = caughey_etal_2017_res$coef[3],
  "se_est" = caughey_etal_2017_res$se[3],
  "p_val" = caughey_etal_2017_res$pv[3]
)

#(2) carson_sievert_2017: 

carson_sievert_2017_df <- read.csv("data/cl-replication-data/carson_sievert_2017.csv")
carson_sievert_2017_df <- carson_sievert_2017_df[carson_sievert_2017_df$margin < 5, ]

carson_sievert_2017_fit <- lm(
  dpres ~ I(rv)*dem_seat + I(rv^2)*dem_seat,
  data = carson_sievert_2017_df
)

carson_sievert_2017_vcov <- sandwich::vcovHC(carson_sievert_2017_fit, type="HC2")

rep_res <- rbind(rep_res, data.frame(
  "estimate_code" = "carson_sievert_2017",
  "beta_est" = carson_sievert_2017_fit$coefficients["dem_seat"],
  "se_est" = sqrt(carson_sievert_2017_vcov["dem_seat", "dem_seat"]),
  "p_val" = NA
))

#(3) dbk_warshaw_2016:

dbk_warshaw_2016_df <- read.csv("data/cl-replication-data/dbk_warshaw_2016.csv")

dbk_warshaw_2016_bw <- rdbwselect_2014(
  y = dbk_warshaw_2016_df$Total.Expenditure.D2,
  x = dbk_warshaw_2016_df$demshare,
  kernel = "uni"
) 

dbk_warshaw_2016_res <- rdrobust(
  y = dbk_warshaw_2016_df$Total.Expenditure.D2 , 
  x = dbk_warshaw_2016_df$demshare, 
  h = dbk_warshaw_2016_bw$bws[1], 
  b = dbk_warshaw_2016_bw$bws[2], 
  all = TRUE, 
  kernel = "uni"
)

rep_res <- rbind(rep_res, data.frame(
  "estimate_code" = "dbk_warshaw_2016",
  "beta_est" = dbk_warshaw_2016_res$coef[1],
  "se_est" = dbk_warshaw_2016_res$se[1],
  "p_val" = dbk_warshaw_2016_res$pv[1]
))

# Save cleaned CSV file
write.csv(rep_res, "data/meta-replication.csv", row.names = FALSE)
